Para evaluar el impacto de las observaciones asimiladas se van a comparar dos experimentos de asimilación:
Las siguientes figuras muestran la diferencia en la humedad absoluta y la temperatura entre los dos experimentos a lo largo de todo el periodo de asimilación. Se realizó una interpolación a niveles de presión y luego se calculó un promedio espacial para obtener un perfil vertical en cada tiempo.
La diferencia en la humedad absoluta es importante solo en niveles bajos por debajo de 850-900 hPa. Esto parece ser razonable ya que la mayoría de las observaciones asimiladas son de superficie. En general CONV+AUT tiene valores de humedad más altos cerca de superficie. Por el contrario, la diferencia en la temperatura es importante en todos los niveles. En niveles bajos la diferencia tiene una especie de ciclo diario que se suaviza cuando la convección empieza a desarrollarse en el dominio. También llama la atención el calentamiento en niveles medios/altos durante el 22 y un enfriamientos por arriba.
Será posible que la convección esté ayudando a propagar el efecto de las observaciones en superficie? Porque no se ve en la humedad?
Otra manera de analizar el impacto que generan las observaciones es calculando el “update” o sea la diferencia entre el análisis y el guess (que en este caso es un pronóstico a una hora).
Y acá aparecen cosas raras. PORQUE EL UPDATE DE EXP CONV ES CASI CERO!!! Entre las 18 de 20/11 y las ~04 del 21/11 hay una diferencia del orden de 1e-4, después de eso es CERO y ocurre tanto para la humedad com para la temperatura. 😭😭😭😭😭
Voy a investigar donde metí la pata esta vez.
Si observamos la distribución espacial de la humedad en distintos tiempos (promediando los niveles de presión entre 1000 y 900 hPa), vemos que en general el experimento CONV+AUT genera un entorno más húmedo hasta el sur de Córdoba y centro de Buenos Aires.
Será porque el viento del norte es más intenso en este experimento?
Pero si observamos el campo de la diferencia del viento en el nivel más bajo en la siguiente figura, en general, el viento del norte es menos intenso en CONV+AUT (se observa como viento del sur en la figura pero en realidad indica que la magnitud del viento en CONV+AUT es menor que en CONV). De hecho a las 06Z del 22/11 el viento en niveles bajos pierde intensidad en la región central del dominio.
En las tres variables analizadas se observa una zona al sur de Córdoba y Buenos Aires (donde podría estar ubicado el frente, habría que revisar esto) donde las diferencias entre los dos experimentos es muy importante. Si comparamos estas figuras con la imagen de IR del GOES-16, vemos que la linea donde se genera la convección se ubica donde las diferencias entre los análisis es mayor. Si bien solo con esto es imposible saber si el impacto de las observaciones de superficie fue positivo (o sea que acerca la simulación al estado real de la atmósfera), es interesante ver que el impacto se da justo donde ocurren los procesos que nos interesan.
Es posible que en los experimentos la convección se desarrolle en lugares ligeramente distintos o desfazados en el tiempo y eso genera las diferencias que se ven el gráfico. Sería interesante intentar ver donde hay convección en cada caso
# map <- rnaturalearth::ne_states(country = c("argentina", "Brazil", "Chile", "Uruguay", "Paraguay", "Bolivia"), returnclass = "sf")
map <- sf::read_sf("mapa/mapa.shp")
geom_mapa <- function() {
geom_sf(data = map, fill = NA, color = "black", size = 0.2, inherit.aes = FALSE)
}
dbz <- ReadNetCDF("analisis/dbz/maxdbz_ana_E5_20181122060000.nc", vars = c(dbz = "max_dbz", lon = "XLONG", lat = "XLAT")) %>%
.[, exp := "CONV+AUT"]
dbz_E4 <- ReadNetCDF("analisis/dbz/maxdbz_ana_E4_20181122060000.nc", vars = c(dbz = "max_dbz", lon = "XLONG", lat = "XLAT")) %>%
.[, exp := "CONV"]
pal <- wesanderson::wes_palette("Zissou1", 100, type = "continuous")
rbind(dbz, dbz_E4) %>%
.[dbz > 15] %>%
ggplot(aes(lon, lat)) +
geom_point(aes(color = dbz), size = 0.3) +
scale_color_gradientn(colors = pal) +
geom_mapa() +
coord_sf(xlim = range(dbz$lon), ylim = range(dbz$lat)) +
labs(title = "Reflectividad máxima mayor a 15 dbz",
subtitle = "2018-11-22 06:00:00") +
facet_wrap(~exp) +
theme_minimal()
Vamos a reconocer que hay un elefante en la sala y que algo en el sistema de asimilación no está funcionando como debería. Vamos a analizar las observaciones asimiladas por cada experimento.
| exp | usage.flag | ps | q | t | uv |
|---|---|---|---|---|---|
| E5 | -1 | 230444 | 1318 | 57160 | 164318 |
| E5 | 1 | 46762 | 64974 | 254253 | 154822 |
De esta figura me preocupan varias cosas:
Respecto del segundo punto: Todas las observaciones tipo 181 de humedad son rechazadas. Resulta que me olvide de sacar el -1 en el archivo convinfo. Parece que hice lo mismo para uv (tipo 287).
Muchas estaciones tienen datos faltantes, no siempre pero en algunos casos llegan a sumar 500 observaciones durante el experimento. Calculando la cantidad de observaciones por ciclo de asimilación da un promedio de 1800, aproximadamente lo que da en el gráfico.
Los prepbufrs contienen 3 tipos de observaciones provenientes de aviones:
De los 3 tipos, 131/231 tiene más observaciones.
## Using 'N' as value column. Use 'value.var' to override
## type -1 1
## 1: 130 1 4
## 2: 131 27461 3793
## 3: 133 62 323
## 4: 230 1 3
## 5: 231 27475 3772
## 6: 233 36 348
Por ahora según la configuración de la errtable algunas observaciones de humedad no se asimilan (también está así en el archivo convinfo). Pero serán tan malas?
| tipo | temp | q | uv |
|---|---|---|---|
| 130 | 2.5 | NA | - |
| 131 | 1.47 | NA | - |
| 133 | 1.47 | 1.94 | - |
| 230 | - | - | 3.6 |
| 231 | - | - | 3.0 |
| 233 | - | - | 2.5 |
Veamos que pinta tiene la diferencia entre las observaciones y el guess.
## Warning: Groups with fewer than two data points have been dropped.
obs[type %in% c("131", "231") & exp == "E5"] %>%
ggplot(aes(obs.guess, pressure)) +
geom_point(alpha = 0.01) +
scale_y_level() +
facet_grid(type ~ usage.flag)
obs[type == 131 & exp == "E5" & date == "2018-11-22 00:00:00" & usage.flag == 1]
## var stationID type dhr lat lon pressure usage.flag flag.prep obs
## 1: t AG1002 131 -0.02 -31.85 297.20 250 1 0 225
## 2: t AG1002 131 -0.13 -32.30 298.05 250 1 0 225
## 3: t AG1000 131 0.00 -39.45 299.08 227 1 0 221
## 4: t AG1000 131 -0.12 -38.58 299.42 227 1 0 222
## 5: t AG1000 131 -0.23 -37.75 299.82 228 1 0 222
## 6: t AG1000 131 -0.52 -35.88 301.22 273 1 0 231
## 7: t AG1000 131 -0.52 -35.82 301.27 282 1 0 232
## 8: t AG1000 131 -0.35 -36.98 300.40 228 1 0 220
## 9: t AG1000 131 -0.43 -36.35 300.87 236 1 0 222
## 10: t AG1000 131 -0.47 -36.13 301.03 245 1 0 225
## 11: t AG1000 131 -0.48 -36.03 301.10 254 1 0 226
## 12: t AG1002 131 -0.52 -33.97 300.77 401 1 0 248
## 13: t AG1002 131 -0.55 -34.10 300.98 472 1 0 259
## 14: t AG1002 131 -0.57 -34.17 301.13 532 1 0 267
## 15: t AG1002 131 -0.57 -34.13 301.05 500 1 0 263
## 16: t AG1002 131 -0.58 -34.22 301.20 563 1 0 271
## 17: t AG1002 131 -0.50 -33.92 300.68 386 1 0 245
## 18: t AG1002 131 -0.50 -33.87 300.62 370 1 0 243
## 19: t AG1002 131 -0.40 -33.43 299.98 277 1 0 232
## 20: t AG1002 131 -0.42 -33.53 300.15 304 1 0 237
## 21: t AG1002 131 -0.42 -33.48 300.07 286 1 0 234
## 22: t AG1002 131 -0.43 -33.60 300.23 313 1 0 238
## 23: t AG1002 131 -0.45 -33.72 300.40 338 1 0 242
## 24: t AG1002 131 -0.45 -33.67 300.32 326 1 0 240
## 25: t AG1002 131 -0.47 -33.77 300.47 349 1 0 242
## 26: t AG1002 131 -0.48 -33.82 300.55 358 1 0 243
## 27: t AG1002 131 -0.50 -33.92 300.68 386 1 0 245
## 28: t AG1002 131 -0.50 -33.87 300.62 370 1 0 243
## 29: t AG1002 131 -0.25 -32.73 298.93 250 1 0 226
## 30: t AG1000 131 -0.53 -35.75 301.32 295 1 0 235
## 31: t AG1000 131 -0.55 -35.68 301.37 308 1 0 237
## 32: t AG1000 131 -0.57 -35.60 301.42 319 1 0 240
## 33: t AG1000 131 -0.57 -35.53 301.47 338 1 0 243
## 34: t AG1000 131 -0.60 -35.40 301.57 368 1 0 246
## 35: t AG1000 131 -0.60 -35.33 301.62 384 1 0 248
## 36: t AG1000 131 -0.62 -35.25 301.67 405 1 0 251
## 37: t AG1000 131 -0.63 -35.18 301.72 422 1 0 252
## 38: t AG1000 131 -0.63 -35.12 301.73 443 1 0 254
## 39: t AG1000 131 -0.65 -35.03 301.77 466 1 0 257
## 40: t AG1000 131 -0.67 -34.97 301.80 494 1 0 261
## 41: t AG1000 131 -0.67 -34.90 301.83 523 1 0 266
## 42: t AG1000 131 -0.68 -34.83 301.87 556 1 0 270
## 43: t AG1000 131 -0.70 -34.77 301.90 579 1 0 271
## 44: t AG1000 131 -0.70 -34.70 301.90 600 1 0 274
## 45: t AG1000 131 -0.72 -34.65 301.85 639 1 0 278
## 46: t AG1000 131 -0.73 -34.60 301.82 691 1 0 284
## 47: t AG1000 131 -0.73 -34.55 301.77 747 1 0 289
## 48: t AG1000 131 -0.77 -34.50 301.60 925 1 0 297
## 49: t AG1000 131 -0.77 -34.50 301.65 862 1 0 292
## 50: t AG1000 131 -0.78 -34.53 301.57 963 1 0 300
## 51: t AG1000 131 -0.78 -34.52 301.57 949 1 0 299
## 52: t AG1000 131 -0.80 -34.55 301.57 981 1 0 301
## 53: t AG1002 131 -0.60 -34.28 301.33 637 1 0 278
## 54: t AG1002 131 -0.60 -34.25 301.27 600 1 0 275
## 55: t AG1002 131 -0.62 -34.32 301.40 673 1 0 281
## 56: t AG1002 131 -0.63 -34.37 301.52 744 1 0 288
## 57: t AG1002 131 -0.63 -34.35 301.47 691 1 0 284
## 58: t AG1002 131 -0.65 -34.40 301.57 782 1 0 292
## 59: t AG1002 131 -0.67 -34.43 301.58 865 1 0 292
## 60: t AG1002 131 -0.70 -34.52 301.55 984 1 0 300
## var stationID type dhr lat lon pressure usage.flag flag.prep obs
## obs.guess obs2 obs.guess2 rerr exp date
## 1: -1.8900 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 2: -1.0800 1e+10 187 1.00 E5 2018-11-22 03:00:00
## 3: 1.5800 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 4: 1.5000 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 5: 1.2700 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 6: 1.3100 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 7: 1.2200 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 8: -0.6940 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 9: -0.2530 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 10: 0.2790 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 11: 0.4190 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 12: -2.0900 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 13: 0.1620 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 14: 1.5900 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 15: 0.9630 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 16: 2.4800 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 17: -2.9700 1e+10 131 1.41 E5 2018-11-22 03:00:00
## 18: -2.2400 1e+10 131 1.41 E5 2018-11-22 03:00:00
## 19: 0.7330 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 20: 1.4400 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 21: 1.3900 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 22: 0.9950 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 23: 0.9480 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 24: 1.0600 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 25: -0.5290 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 26: -1.3000 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 27: -2.9700 1e+10 131 1.41 E5 2018-11-22 03:00:00
## 28: -2.2400 1e+10 131 1.41 E5 2018-11-22 03:00:00
## 29: -1.0700 1e+10 187 1.00 E5 2018-11-22 03:00:00
## 30: 1.4500 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 31: 1.3000 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 32: 1.4200 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 33: 1.1500 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 34: -0.4040 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 35: -0.8590 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 36: -1.2800 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 37: -2.3100 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 38: -2.5400 1e+10 181 1.00 E5 2018-11-22 03:00:00
## 39: -0.9810 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 40: 0.3660 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 41: 1.7200 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 42: 2.2800 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 43: 1.4300 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 44: 1.8800 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 45: 1.0400 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 46: 1.1500 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 47: 2.0100 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 48: -0.1180 1e+10 131 1.29 E5 2018-11-22 03:00:00
## 49: -0.2900 1e+10 131 1.15 E5 2018-11-22 03:00:00
## 50: 0.2160 1e+10 131 1.38 E5 2018-11-22 03:00:00
## 51: -0.0554 1e+10 131 1.35 E5 2018-11-22 03:00:00
## 52: 0.3250 1e+10 131 1.43 E5 2018-11-22 03:00:00
## 53: 0.9660 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 54: 1.7300 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 55: 0.5010 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 56: 1.2800 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 57: 1.4300 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 58: 3.6400 1e+10 131 1.00 E5 2018-11-22 03:00:00
## 59: -0.0714 1e+10 131 1.15 E5 2018-11-22 03:00:00
## 60: -0.4960 1e+10 181 1.43 E5 2018-11-22 03:00:00
## obs.guess obs2 obs.guess2 rerr exp date